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ABSTRACT 

Acceleration  and  transport  of  high-energy  particles  and  fluid  dynamics  of  atmospheric  plasma  are  interrelated  aspects 
of  solar  flares,  but  for  convenience  and  simplicity  they  were  artificially  separated  in  the  past.  We  present  here  self- 
consistently  combined  Fokker-Planck  modeling  of  particles  and  hydrodynamic  simulation  of  flare  plasma.  Energetic 
electrons  are  modeled  with  the  Stanford  unified  code  of  acceleration,  transport,  and  radiation,  while  plasma  is 
modeled  with  the  Naval  Research  Laboratory  flux  tube  code.  We  calculated  the  collisional  heating  rate  directly  from 
the  particle  transport  code,  which  is  more  accurate  than  those  in  previous  studies  based  on  approximate  analytical 
solutions.  We  repeated  the  simulation  of  Mariska  et  al.  with  an  injection  of  power  law,  downward-beamed  electrons 
using  the  new  heating  rate.  For  this  case,  a  ~10%  difference  was  found  from  their  old  result.  We  also  used  a  more 
realistic  spectrum  of  injected  electrons  provided  by  the  stochastic  acceleration  model,  which  has  a  smooth  transition 
from  a  quasi-thermal  background  at  low  energies  to  a  nonthermal  tail  at  high  energies.  The  inclusion  of  low-energy 
electrons  results  in  relatively  more  heating  in  the  corona  (versus  chromosphere)  and  thus  a  larger  downward  heat 
conduction  flux.  The  interplay  of  electron  heating,  conduction,  and  radiative  loss  leads  to  stronger  chromospheric 
evaporation  than  obtained  in  previous  studies,  which  had  a  deficit  in  low-energy  electrons  due  to  an  arbitrarily 
assumed  low-energy  cutoff.  The  energy  and  spatial  distributions  of  energetic  electrons  and  bremsstrahlung  photons 
bear  signatures  of  the  changing  density  distribution  caused  by  chromospheric  evaporation.  In  particular,  the  density 
jump  at  the  evaporation  front  gives  rise  to  enhanced  emission,  which,  in  principle,  can  be  imaged  by  X-ray  telescopes. 

This  model  can  be  applied  to  investigate  a  variety  of  high-energy  processes  in  solar,  space,  and  astrophysical  plasmas. 
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Sun:  X-rays,  gamma  rays 
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1.  INTRODUCTION 

A  solar  flare,  as  one  of  the  most  prominent  manifestations  of 
solar  activity,  has  many  faces  among  which  are  acceleration  and 
transport  of  high-energy  particles  and  the  dynamic  response 
of  atmospheric  plasma.  It  is  generally  believed  that  magnetic 
reconnection  in  the  corona  is  the  primary  energy  release  mech¬ 
anism  that  leads  to  plasma  heating  and  particle  acceleration. 
The  heated  plasma  and  accelerated  particles  (primarily  elec¬ 
trons)  produce  bremsstrahlung  X-rays  at  the  apex  of  the  flare 
loop  observed  as  a  loop-top  (LT)  source  (e.g.,  Masuda  et  al. 
1994;  Petrosian  et  al.  2002;  Krucker  &  Lin  2008).  Some  of 
the  released  energy  is  transported  down  the  closed  magnetic 
loop  by  nonthermal  particles  (electrons  and  ions)  and  thermal 
conduction,  which  contribute  to  energy  gain  in  various  layers 
of  the  atmosphere.  Electrons  give  up  most  of  their  energy  to 
ambient  particles  via  Coulomb  collisions,  and  produce  hard 
X-rays  (HXRs)  primarily  at  the  footpoints  (FPs)  of  the  loop  in 
the  dense  transition  region  (TR)  and  chromosphere  (see  Hoyng 
et  al.  1981;  Sakao  1994;  Saint-Hilaire  et  al.  2008).  Accelerated 
protons  and  heavier  ions,  on  the  other  hand,  cause  nuclear  re¬ 
actions  while  colliding  with  background  particles  and  produce 
y-rays.  Some  accelerated  electrons  and  ions  escape  along  open 
magnetic  field  lines  into  interplanetary  space  and  are  observed 
as  solar  energetic  particles  by  in  situ  instruments  (e.g.,  Lin  1985; 
Liu  et  al.  2004b;  Krucker  et  al.  2007).  The  rate  of  energy  gain 
in  the  chromosphere,  if  exceeding  the  combined  local  radiative 
and  conductive  cooling  rate,  can  rapidly  heat  the  plasma  up  to 
a  temperature  of  ~10^  K.  The  resulting  overpressure  drives  an 


upward  mass  flow  at  a  speed  up  to  hundreds  of  km  s“',  which 
fills  the  flare  loop  with  a  hot  plasma,  giving  rise  to  the  gradual 
increase  of  soft  X-ray  (SXR)  emission.  This  process,  termed 
chromospheric  evaporation  by  Neupert  (1968),  can  influence 
particle  transport  by  changing  the  ambient  density  in  the  loop 
on  timescales  of  tens  of  seconds,  and  affect  heat  conduction 
by  changing  the  loop  temperature  distribution  at  the  same  time. 
Collisional  and  conductive  heating  will  be  consequently  mod¬ 
ified,  and  in  turn,  so  will  the  dynamic  atmospheric  response. 
On  longer  timescales  of  minutes,  as  magnetic  reconnection  pro¬ 
ceeds  and  new  loops  are  formed  or  excited,  the  above  processes 
repeat  sequentially  in  newer  loops. 

1.1.  Motivation  for  This  Study 

The  aforementioned  processes  are  coupled  in  a  circular  chain, 
but  due  to  great  complexity  of  the  subject,  previous  researchers 
tended  to  focus  on  one  process  at  a  time  while  assuming  some 
simple  forms  for  others.  Past  investigations  fall  into  two  general 
categories:  (1)  acceleration  and  transport  of  particles  and  (2) 
fluid  dynamics  of  atmospheric  plasma. 

Various  mechanisms  have  been  proposed  for  particle  ac¬ 
celeration.  Among  the  agents  of  acceleration  are  DC  elec¬ 
tric  fields  (Holman  1985;  Litvinenko  1996;  Zharkova  & 
Gordovskyy  2004),  shocks  (Tsuneta  &  Naito  1998),  and  tur¬ 
bulence  or  plasma  waves  (Ramaty  1979;  Hamilton  &  Petrosian 
1992;  Miller  et  al.  1996;  Petrosian  &  Liu  2004).  Particle  trans¬ 
port  is  comparatively  better  understood  and  previous  studies 
usually  assumed  a  hydrostatic  atmosphere.  Early  analytical 
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studies  (Brown  1973;  Petrosian  1973;  Lin  &  Hudson  1976; 
Emslie  1978;  Chandrashekar  &  Emslie  1987)  took  various  ap¬ 
proximations  (e.g.,  neglecting  pitch-angle  diffusion)  to  allow 
the  problem  to  be  tractable.  The  numerical  study  of  Leach 
&  Petrosian  (1981)  improved  on  this  by  solving  the  Fokker- 
Planck  transport  equation  with  inclusion  of  pitch-angle  changes 
due  to  Coulomb  collision  and  magnetic  mirroring.  This  was 
later  extended  to  the  relativistic  regime  by  including  energy 
losses  and  pitch-angle  changes  due  to  synchrotron  radiation 
(McTieman  &  Petrosian  1990).  Similar  Fokker-Planck  stud¬ 
ies  of  particle  transport  were  performed  by  MacKinnon  & 
Craig  (1991),  McClements  (1992),  and  Syniavskii  &  Zharkova 
(1994). 

Fluid  dynamics  of  the  magnetized  atmosphere  in  response 
to  flare  heating  can  be  best  studied  with  a  multi-dimensional 
magnetohydrodynamic  (MHD)  model,  but  for  simplicity  most 
efforts  were  invested  in  one-dimensional  hydrodynamic  (HD) 
simulations.  This  is  justified  for  the  solar  corona  where  the 
magnetic  pressure  dominates  the  gas  pressure  (low-/!  plasma) 
and  resistivity  is  low.  As  a  result,  plasma  is  only  allowed  to  flow 
along  the  magnetic  field  lines  due  to  the  line-tying  condition. 
Previous  HD  models  (MacNeice  et  al.  1984;  Nagai  &  Emslie 
1984;  Emslie  &  Nagai  1985;  Fisher  et  al.  1985a,  1985b,  1985c; 
Mariska  et  al.  1989;  Gan  &  Fang  1990;  Emslie  et  al.  1992) 
usually  assumed  a  power-law  spectrum  of  accelerated  electrons 
injected  at  the  apex  of  the  loop,  and  calculated  collisional  heating 
along  the  loop  by  these  electrons  from  approximate  analytical 
solutions  of  particle  transport  mentioned  above.  Abbett  & 
Hawley  (1999)  and  Allred  et  al.  (2005)  improved  on  previous 
studies  by  including  detailed  calculation  of  radiative  transfer  in 
the  atmosphere. 

There  are  theoretical  and  observational  motivations  to  inves¬ 
tigate  the  particle  and  fluid  aspects  of  a  solar  flare  together  in  a 
self-consistent  manner.  From  a  theoretical  point  of  view,  such  an 
investigation  is  demanded  in  order  to  retrieve  missing  physics 
when  the  two  aspects  were  studied  separately.  It  has  also  become 
technically  more  feasible,  thanks  to  advances  in  both  aspects 
over  the  last  three  decades  and  particularly  in  recent  years.  Sev¬ 
eral  independent  studies  (Miller  &  Mariska  2005;  Winter  et  al. 
2007;  Winter  2009)  along  this  direction  are  already  under  way, 
but  none  of  them  has  been  completed.  From  an  observational 
point  of  view,  new  observations,  particularly  X-ray  images  and 
spectra  obtained  by  the  current  Ramaty  High  Energy  Solar  Spec¬ 
troscopic  Imager  {RHESST)  and  previous  Yohkoh  missions,  have 
posed  challenges  to  the  existing  theories.  For  example,  in  recent 
studies  of  the  Neupert  (1968)  effect,  Veronig  et  al.  (2005)  and 
Liu  et  al.  (2006a)  found  that,  unexpectedly,  the  nonthermal  elec¬ 
tron  energy  deposition  power,  which  is  more  physically  related 
to  the  plasma  thermal  energy  change  rate,  did  not  yield  a  better 
correlation  with  the  time  derivative  of  the  SXR  flux  than  the  con¬ 
ventional  HXR  flux  did.  In  an  event  of  chromospheric  evapora¬ 
tion  imaged  by  RHESSI  for  the  first  time,  Liu  et  al.  (2006a)  found 
X-ray  sources  moving  from  the  FPs  to  the  LT  at  very  high  speeds 
(~10^  km  s“^).  More  interestingly,  Sui  et  al.  (2006)  found  dou¬ 
ble  nonthermal  sources  moving  first  downward  from  the  LT 
toward  the  FPs  and  then  upward  along  the  loop.  To  fully  under¬ 
stand  these  observations  requires  a  joint  study  of  acceleration 
and  transport  of  particles  and  fluid  dynamics  of  the  atmospheric 
response. 

1.2.  Approach  of  This  Study 

With  the  goal  to  investigate  the  coupled  processes  of  acceler¬ 
ation,  transport,  and  hydrodynamics  in  solar  flares,  we  present 


here  combined  Fokker-Planck  modeling  of  particles  and  HD 
simulation  of  plasma.  (1)  The  Fokker-Planck  model  utilizes 
the  Stanford  unified  code  of  particle  acceleration,  transport, 
and  bremsstrahlung  radiation  (Petrosian  et  al.  2001).  The  trans¬ 
port  and  radiation  calculation  is  based  on  the  work  of  Leach  & 
Petrosian  (1981,  1983)  and  McTiernan  &  Petrosian  (1990).  The 
acceleration  module  of  the  code  adopts  the  stochastic  accelera¬ 
tion  model  of  Petrosian  &  Liu  (2004,  hereafter  PL04),  which  has 
inherited  knowledge  accumulated  over  a  decade  (Hamilton  & 
Petrosian  1992;  Dung  &  Petrosian  1994;  Park  &  Petrosian  1995, 
1996;  Park  et  al.  1997).  When  compared  with  observations,  this 
model  has  many  attractive  features  and  advantages  (Liu  et  al. 
2004b,  2006b,  2008)  over  other  mechanisms.  (2)  The  HD  sim¬ 
ulation  uses  the  Naval  Research  Laboratory  (NRL)  Solar  Flux 
Tube  Model  (Mariska  et  al.  1989,  hereafter  MEL89),  which,  as 
a  modified  version  of  the  Mariska  et  al.  (1982)  model,  provides 
excellent  treatment  of  fluid  dynamics  and  has  been  widely  used 
in  studying  atmospheric  response  to  flare  heating  (e.g.,  Warren 
&  Doschek  2005). 

One  of  the  major  advances  marked  by  this  study  is  the  more 
accurate  and  self-consistent  evaluation  of  the  collisional  heating 
rate  by  nonthermal  electrons.  This  heating  rate  is  critical  to  HD 
simulation  of  flares,  but  was  not  properly  calculated  previously 
in  two  major  aspects.  (1)  The  calculation  of  energy  loss  of 
energetic  electrons  and  thus  the  heating  rate  was  based  on 
approximate  analytical  solutions  (e.g.,  Brown  1973;  Emslie 
1978),  which  incorporated  only  pitch-angle  growth  due  to 
Coulomb  collisions,  but  in  reality  the  pitch-angle  change  is 
a  diffusion  process.  This  will  be  remedied  in  this  study  with 
the  inclusion  of  a  full  Fokker-Planck  treatment  of  electron 
transport.  (2)  Another  previous  drawback  was  the  use  of  an 
unrealistic  spectrum  of  injected  electrons,  which  usually  was  a 
power  law  with  a  low-energy  cutoff.  Fisher  et  al.  (1985c),  for 
example,  assumed  a  sharp  low-energy  cutoff  at  E^  =  20  keV 
(i.e.,  no  electrons  below  £c),  while  MEL89  introduced  a  “soft” 
cutoff  below  which  the  spectrum  is  still  a  power  law  with  a 
positive  slope.  It  should  be  noted  that  rather  than  an  intrinsic 
property  of  the  primary  accelerated  electron  population,  a  low- 
energy  cutoff  or  turnover  in  the  electron  spectrum  inferred 
from  X-ray  observations  can  result  from  secondary  effects 
such  as  return  currents  (Zharkova  &  Gordovskyy  2006)  and 
photospheric  albedo  (Danger  &  Petrosian  1977;  Bai  &  Ramaty 
1978;  Sui  et  al.  2007).  The  collisional  heating  rate  is  sensitive 
to  the  injected  electron  spectrum  and  thus  the  use  of  an 
incorrect  spectrum  would  make  the  HD  simulation  deviate 
from  reality  significantly.  PL04  has  provided  a  more  realistic 
electron  spectrum  that  has  a  continuous  span  from  a  quasi- 
thermal  distribution  at  low  energies  to  a  nonthermal  tail  at  high 
energies,  avoiding  an  unnecessary  low-energy  cutoff.  It  also 
gives  good  fits  to  both  LT  and  FP  X-ray  spectra  obtained  by 
RHESSI.  Such  an  electron  spectrum  is  used  in  this  work.  As 
we  will  see  later,  the  low-energy  electrons,  which  otherwise 
would  have  been  missing  if  a  cutoff  were  to  be  present,  play  an 
important  role  in  heating  and  in  influencing  the  subsequent  HD 
evolution. 

We  present  the  numerical  model  in  Section  2,  techniques 
to  combine  different  modules  of  the  model  in  Section  3,  and 
simulation  runs  in  Section  4.  We  compare  the  HD  characteris¬ 
tics  of  different  simulations  in  Section  5  and  examine  the  HD 
effects  on  particle  transport  and  X-ray  emission  in  Section  6. 
Conclusions  and  discussions  are  given  in  Section  7.  This  model 
is  based  on  the  PhD  thesis  of  Liu  (2006),  whose  revised  edi¬ 
tion  has  also  appeared  as  a  book  (Liu  2008).  The  model  has 
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Figure  1.  Top:  geometry  of  the  model  flare  loop.  Q  is  the  pitch  angle  of  the 
electron  with  a  velocity  Vg  in  the  guiding  magnetic  field  Bq.  Bottom:  initial 
distribution  of  logarithmic  temperature  T,  gas  pressure  P  (left  scale),  and  electron 
number  density  rie  (right  scale)  vs.  distance  in  the  transport  region  of  the  loop. 
Pressure  is  scaled  upward  by  a  factor  of  100. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


been  refined  ever  since  and  the  numerical  results  presented  here 
are  new.  Our  new  calculations  and  analyses  include:  (1)  re¬ 
computing  the  MEL89  simulation  using  our  transport  code  (see 
Section  4.2);  (2)  detailed  analysis  of  the  interplay  of  heating  and 
cooling  (Section  5.2);  and  (3)  examining  the  temperature  distri¬ 
bution  of  plasma  velocity  which  can  directly  be  compared  with 
Doppler  observations  (Section  5.3).  In  order  to  achieve  more  ac¬ 
curacy,  we  have  extended  the  number  of  pitch-angle  bins  from  24 
to  100. 


the  loop.  Parts  (1),  (2),  and  (4)  are  inherited  from  the  Stanford 
unified  particle  code  (Petrosian  et  al.  2001)  in  which  the  accel¬ 
eration  module  has  been  revised  according  to  PL04,  while  part 
(3)  is  adopted  from  the  NRL  flux  tube  code  (MEL89).  Details 
of  the  four  parts  are  described  in  the  following  subsections. 

Note  that  sequential  energizing  (or  excitation)  propagating 
from  one  flare  loop  to  another  has  been  observed  as  the  apparent 
motions  of  HXR  LT  (Gallagher  et  al.  2002;  Sui  &  Holman 
2003)  and  FP  (Grigis  &  Benz  2005;  Yang  et  al.  2009;  Liu 
et  al.  2009)  sources.  In  our  simulations,  we  set  the  duration 
of  the  impulsive  phase  to  be  60  s.  This  time,  according  to 
Schrijver  et  al.  (2006),  translates  into  an  apparent  FP  speed  of 
^  =  10  km  s“'  (where  r  —  0.3  Mm),  which  is  comparable 
to  the  observed  values  (Liu  et  al.  2004a).  Our  single  loop 
scenario  is  thus  legitimate  within  the  relevant  timescale  and 
can  be  viewed  as  an  elementary  process  of  sequential  excitation 
of  multiple  loops.  Evolution  on  longer  timescales  (say,  >100  s) 
involves  multiple  loops  and  can  be  studied  by  superimposing 
sequential  single-loop  simulations  (e.g.,  Warren  2006). 


2.7.  Stochastic  Acceleration 

The  stochastic  acceleration  model  of  PL04  addresses  electron 
and  proton  acceleration  by  plasma  waves  propagating  parallel 
to  the  background  magnetic  field  Bq.  According  to  this  model, 
large-scale  turbulence  or  long-wavelength  plasma  waves  are 
generated  in  the  corona  as  a  result  of  magnetic  reconnection. 
The  turbulence,  cascading  to  smaller  scales,  heats  plasma  and 
accelerates  particles  in  a  region  near  the  top  of  the  flaring 
loop.  The  heated  plasma  and  accelerated  electrons  produce  the 
observed  thermal  and  nonthermal  X-rays,  respectively,  in  the 
acceleration  region  or  the  LT  source  (Xu  et  al.  2008;  Liu  et  al. 
2008).  Here  we  briefly  repeat  the  mathematical  description  of 
the  model.  The  Fokker-Planck  equation  that  governs  electron 
acceleration,  in  general,  can  be  written  as 


dt 
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2.  SIMULATION  MODEL 

Here,  we  consider  the  dynamic  evolution  of  a  single  flare  loop 
perpendicular  to  the  solar  surface.  As  shown  in  Figure  1  (top), 
an  acceleration  region  of  length  L  =  5  Mm  is  located  at  the  top 
of  the  model  loop,  sandwiched  between  two  symmetric  quarter- 
circles  called  the  transport  region  of  length  Smax  =  14  Mm. 
The  loop  has  a  uniform  circular  cross  section  a(s)  —  ao  —  itr^ 
with  a  constant  radius  of  r  =  0.3  Mm  at  any  distance  s  from 
the  edge  of  the  acceleration  region.  In  its  initial  state  (Figure  1 , 
bottom)  the  loop  spans  from  the  hot  {T  >  10®  K),  tenuous 
corona  to  the  cold  (T  =  10"^  K),  dense  chromosphere,  with  the 
TR  (defined  here  as  the  lowest  point  where  r  ^  10®  K)  located 
at  ^tr  —  10  Mm. 

The  simulation  model  includes  four  parts:  (1)  the  stochas¬ 
tic  acceleration  code  generates  a  spatially  averaged  distribution 
(in  energy)  of  high-energy  electrons  in  the  acceleration  region. 
This  distribution  is  fed  to  the  transport  region  where  (2)  the 
transport  code  computes  the  electron  distribution  (in  energy 
and  pitch  angle)  and  collisional  heating  rate  as  a  function  of 
distance  (depth),  and  (3)  the  hydrodynamic  code  simulates  the 
atmospheric  response  to  this  heating.  Finally,  (4)  the  radiation 
code  calculates  the  corresponding  bremsstrahlung  emission  in 


where  /ac  =  /ac(L  E),  in  units  of  electrons  cm“®  keV“\  is 
the  angle-integrated  and  spatially  averaged  electron  distribution 
function  in  energy  space  (subscript  “ac”  denotes  the  acceleration 
region),  E  is  the  electron  kinetic  energy,  D(E)  is  the  energy 
diffusion  coefficient,  A(E)  the  direct  acceleration  rate,  T^sc{E) 
the  particle  escape  time,  Q(E)  is  the  rate  of  background 
electrons  being  supplied  to  the  acceleration  region  that  serves 
as  a  source  term,  and  finally  Ei  —  EcouX  +  Esynch  is  the  absolute 
value  of  the  net  energy-loss  rate  that  is  a  sum  of  the  Coulomb 
and  synchrotron  loss  rates. 

In  order  to  solve  Equation  (1)  one  must  first  evaluate  all  the 
terms.  The  energy-loss  rates  Ecoui  and  Esynch  are  well  known 
(see  Equation  (18)  of  PL04)  and  the  source  term  Q{E)  is  to  be 
prescribed  with  a  specific  model  (assumed  to  be  a  thermal  or 
Maxwellian  distribution  here).  The  central  task  is  then  to  ob¬ 
tain  D(E),  A(E),  and  T^sdE).  They  are  determined  through  the 
dispersion  relation  of  the  plasma  waves  (PL04,  Equation  (28)) 
and  the  wave-particle  resonance  condition  (PL04,  Equation  (4)). 
Following  PL04,  we  assume  a  fully  ionized  H  and  ^He  plasma 
with  a  relative  abundance  of  electron/proton/a -particle  = 
1/0.84/0.08,  and  a  broken  power-law  spectrum  of  the  turbu¬ 
lence  given  by  Equation  (29)  in  PL04  with  the  relevant  spectral 
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indices  qi  —  2,  q  —  —l.l  (Kolmogorov  value),  and  qu  —  —4. 
The  characteristic  acceleration  rate  given  by  Equation  (30) 
in  PL04  represents  the  rate  of  wave-particle  interaction  and  de¬ 
pends  on  the  level  of  turbulence. 

Once  all  the  coefficients  have  been  evaluated.  Equation  (1)  is 
solved  numerically  using  the  flux  conservative  finite  difference 
scheme  of  Chang  &  Cooper  (1970)  described  in  Park  &  Pet¬ 
rosian  (1996).  Here  we  assume  a  homogeneous  acceleration  re¬ 
gion  and  obtain  a  steady  state  solution  of  /ac(E)  (i.e.,  d/dt  —  0). 
The  angle-integrated  flux  in  the  acceleration  region  is  — 

VeficiE),  where  Vg  is  the  electron  velocity.  We  then  calculate  the 
flux  of  electrons  that  escape  the  acceleration  region  and  enter 
the  transport  region  of  the  flare  loop, 

/ac(£) 

F,sc(E)  =  .  (2) 

^  escv^  / 


RHESSI  flares  (Kontar  &  Brown  2006;  Kasparova  et  al.  2007). 
The  injected  flux  at  the  top  (s  =  0)  of  each  leg  of  the  loop 
is  then  E{E,  /x,  i')|s=o  =  E^sciE)/  /_j  dfx  which  is  equivalent 
to  a  uniform  distribution  in  the  4jt  solid  angle  integrated  over 
the  2jt  range  of  the  azimuthal  angle  0  assuming  axisymmetry. 
With  the  symmetric  assumption,  the  steady-state  calculation  is 
performed  in  only  one  leg  of  the  loop.  We  impose  a  symmetric 
(or  reflective)  boundary  condition  at  5  =  0,  where  a  particle 
leaving  the  computational  domain  is  reflected  back  with  identi¬ 
cal  energy  but  opposite  pitch-angle  cosine,  mimicking  a  particle 
coming  from  the  other  leg  of  the  loop.  (2)  As  to  the  background 
atmosphere,  we  assume  a  fully  ionized  hydrogen  plasma  whose 
distribution  is  taken  from  the  result  of  the  HD  code  described 
next. 

2.3.  Hydrodynamics 


2.2.  Particle  Transport 

The  flux  FgsciE)  is  then  input  to  the  particle  transport  code 
(Leach  &  Petrosian  1981;  McTieman  &  Petrosian  1990)  which 
calculates  the  electron  distribution  in  energy  and  pitch-angle 
space,  and  its  variation  with  distance  while  the  electrons  spiral 
down  magnetic  field  lines  into  deeper  layers  of  the  atmosphere. 
The  code  numerically  solves  the  fully  relativistic,  steady-state, 
Fokker-Planck  equation  (i.e..  Equation  (1)  in  McTiernan  & 
Petrosian  1990),  which  is  similar  to  Equation  (1)  here  and 
includes  energy  loss  (no  energy  diffusion)  and  pitch-angle 
diffusion  due  to  Coulomb  collision,  and  pitch-angle  changes 
due  to  magnetic  mirroring  and  synchrotron  radiation.  Following 
McTieman  (1989),  we  neglect  return  currents  (Syniavskii  & 
Zharkova  1994;  Zharkova  et  al.  1995). 

The  variable^  to  be  solved  in  the  transport  equation  is  the 
electron  flux  spectrum  E{E,iJ,,s)  as  a  function  of  energy 
E,  cosine  pt  =  cos  6  of  pitch  angle  0,  and  distance  s  from 
the  injection  site  at  the  boundary  of  the  acceleration  region. 
E{E ,  pt,  s)dpL  has  units  of  electrons  s“'  cm“^  keV“*  and  is 
evaluated  as 

F(£', /X,  ^)  =  u^/(£', /X,  j)— ,  (3) 

flO 

where  f  (E,  pt,  s)dpL  is  the  number  density  distribution  function 
in  units  of  electrons  cm“^  keV“^  (compare,  the  angle-integrated 
number  density  fac(E)  in  the  above  acceleration  code),  and 
we  integrate  the  differential  electron  flux  Vef(E,  pt,  s)  over  the 
cross-sectional  area  a{s)  of  the  loop  and  then  normalize  it  by  a 
constant  equivalent  area  aq-  As  noted  earlier,  here  we  assume 
a  constant  a{s)  —  aq  for  simplicity,  which  means  a  uniform 
magnetic  field  along  the  loop  and  thus  no  magnetic  mirroring. 

In  addition  to  the  injected  electron  flux  EgsdE)  from  the  ac¬ 
celeration  code,  the  transport  code  requires  the  knowledge  of 
the  ambient  density  and  abundance  distribution  along  the  loop. 
(1)  Here  we  assume  that  Fesc  is  isotropic  in  pitch  angle,  repre¬ 
senting  the  consequence  of  frequent  scatterings  of  electrons  by 
turbulence  in  the  acceleration  region.  This  assumption  is  consis¬ 
tent  with  the  nearly  isotropic,  rather  than  beamed,  distributions 
inferred  from  center-to-limb  variations  of  HXR  and  y-ray  fluxes 
and  spectral  indices  in  observations  obtained  by  the  Solar  Maxi¬ 
mum  Mission  (McTiernan  &  Petrosian  1991),  and  more  recently 
from  atmospheric  albedo  due  to  Compton  back-scattering  in 

^  In  practice,  the  numerical  code  equivalently  solves  for 

F(E,  /X,  s)/ =  c<Pa(s)/aQ,  where  O  is  defined  in  McTieman  &  Petrosian 

(1990)  and /3  =  ujc. 


Hydrodynamics  in  the  transport  region  is  calculated  with  the 
NRL  solar  flux  tube  code  (MEL89)  based  on  Mariska  et  al. 
(1982).  The  code  assumes  a  two-fluid  plasma  composed  of  elec¬ 
trons  and  ions  that  can  only  move  along  the  magnetic  field  in 
a  flux  tube.  The  user-specified  geometry  of  the  tube  (a  uniform 
quarter-circle  in  our  case)  is  characterized  by  the  tube  cross- 
sectional  area  a{s)  and  the  component  of  the  gravitational  ac¬ 
celeration  along  the  tube.  The  code  solves  the  time-dependent, 
one-dimensional  equations  of  mass,  momentum,  and  energy 
conservation  (see  Equations  (l)-(3)  in  MEL89).  The  indepen¬ 
dent  variables  are  the  mass  density  p{s),  fluid  velocity  v{s),  and 
temperature  T{s)  which  we  assume  to  be  the  same  for  electrons 
and  ions.  Because  of  small  masses  of  electrons,  we  neglect 
the  momentum  loss  of  energetic  electrons  to  the  background 
plasma. 

The  volumetric  heating  rate  in  the  energy  equation  is 

5(y)  =  5,(y)  +  Sa  ,  (4) 


where  Se{s)  represents  heating  by  energetic  electrons,  which  is 
provided  by  the  transport  code  (see  Section  3),  and  Sq  =  8.31  x 
10“^  erg  s“'  cm“^  (MEL89)  represents  uniform  background 
heating,  presumably  caused  by  coronal  heating  in  the  quiet  Sun 
active  region.  The  conductive  flux  Econd  and  heating  rate  5cond 
are 


dT 

Tcond  —  F  --  , 

0S 


*5cond  — 


9  Ucond 
ds 


(5) 


where  k  is  the  thermal  conductivity.  The  radiative  energy-loss 
rate  is 


Trad  —  ri^n p^{T')  , 


(6) 


where  and  np  are  the  electron  and  proton  number  density, 
respectively,  which  are  equal  by  our  assumption  of  fully  ionized 
hydrogen  plasma,  and  0(7)  is  the  optically  thin  radiative  loss 
function  (MEL89)  which  has  its  maximum  at  7  ~  (1-3)  x 
10^  K. 

We  select  an  adaptive  mesh  of  450  grids  that  move  with  time 
to  optimize  spatial  resolution  in  the  dense  chromosphere  and 
near  sharp  jumps  at  the  TR  and  evaporation  front.  This  mesh 
is  also  shared  by  the  transport  and  bremsstrahlung  radiation 
codes  in  our  new  model.  A  reflective  (or  symmetric)  boundary 
condition  is  imposed  at  both  the  upper  (loop  apex)  and  lower 
(deep  in  the  chromosphere)  boundaries  of  the  transport  region 
(see  Figure  1),  such  that  the  system  remains  closed. 


2.4.  Bremsstrahlung  Radiation 

Having  obtained  the  electron  flux  from  the  transport  code  and 
the  background  density  from  the  HD  code,  we  calculate  the  thin- 
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target  bremsstrahlung  radiation  intensity  or  photon  emission 
rate,  l(e,  s),  as  a  function  of  photon  energy  e  and  distance  s. 
I(€,  s)  (photons  s“'  cm“^  keV“^)  is  defined  as 


I(€,s) 


da(e,  E) 
de 


(7) 


where  da{e,  E)/de  is  the  angle-averaged  differential 
bremsstrahlung  cross  section  given  by  Koch  &  Motz  (1959), 
and  Finx  —  J  j  E{E,  fi,  s)dti  is  the  angle-integrated  electron 
flux.  Substituting  Fint  with  the  acceleration  region  flux  E^c  gives 
the  LT  emission  lurie),  while  identifying  fint  =  LEihick  yields 
the  spatially  integrated  thick-target  emission  /thick(f)  (Brown 
1971;  Petrosian  1973).  Here  Eihick  =  fc/tMck  is  the  equivalent 
thick-target  electron  flux,  given  by  the  corresponding  number 
density  (Petrosian  &  Donaghy  1999;  PL04) 


/thick(£’)  = 


1  c°° 

LEl  Je 


F,,,{E')dE' . 


(8) 


The  equivalent  FP  emission  is  the  spatially  averaged  pho¬ 
tons  emitted  below  the  TR  (located  at  s^),  7pp(e)  = 
— itr), where  imax  is  the  distance  at  the  lower 
boundary  of  the  loop.  If  the  corona  is  negligibly  tenuous  and  the 
column  depth  at  j^ax  is  large  enough  to  stop  all  HXR  producing 
electrons  of  interest,  (j^ax  —  Sa)Ivp(€)  approaches  /thick(e)- 
Comparison  between  the  HXRs  modeled  here  and  those 
observed  by  the  RHESSI  satellites  can  serve  as  a  unique 
diagnostic  tool  and  will  be  pursued  in  a  future  publication.  Here 
we  show  an  example  (Figure  2)  of  how  well  observed  LT  and  FP 
fluxes  can  be  fitted  with  the  above  equations  using  a  spectrum 
of  accelerated  electrons  given  by  the  PL04  model. 


3.  COMBINING  PARTICLE  AND  HYDRODYNAMIC 
CODES 

The  main  task  in  this  study  is  to  combine  the  Stanford  particle 
code  and  the  NRL  HD  code.  Here  we  assume  that  particle 
acceleration  acts  as  an  independent  driver  of  the  simulation 
and  is  not  affected  by  the  temperature  or  density  evolution  in 
the  transport  region  of  the  loop.  The  task  is  thus  reduced  to 
make  the  transport  module  of  the  particle  code  and  the  HD  code 
communicate  interactively  in  real  time. 

3.1.  Electron  Heating  Rate 


E  (keV) 


Figure  2.  Photon  fluxes  multiplied  by  energy  squared  (e^)  during  the  main 
HXR  peak  of  the  2002  August  3  Xl.O  flare  and  spectral  fits  using  the  stochastic 
acceleration  model  of  PL04.  (a)  Summed  fluxes  of  all  LT  sources  and  all  FP 
sources  (see  Figure  2.6  of  Liu  2006).  (b)  Sum  (squares)  of  the  LT  and  FP  fluxes 
shown  in  (a)  and  sum  (solid  line)  of  the  corresponding  model  fits.  Overplotted 
are  the  spatially  integrated  spectrum  (plus  signs)  and  the  corresponding  preflare 
background  (dotted  line).  The  legend  lists  parameters  used  in  the  model.  (We 
thank  Siming  Liu  and  Yan-Wei  Jiang  for  help  in  producing  this  figure.) 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


and  differentiate  it  to  obtain  the  net  energy  gain  to  the  back¬ 
ground  plasma  in  a  unit  volume. 


The  rate  of  collisional  heating  to  the  background  plasma,  Se 
(erg  s“ '  cm“^),  equals  the  rate  of  energy  loss  from  the  energetic 
electrons.  This  can  be  calculated  from  the  electron  distribution 
given  by  the  transport  code  in  two  equivalent  ways,  the  second 
of  which  is  used  in  this  study. 

Se  can  be  evaluated  from  the  energy-loss  rate  Ecoui  due  to 
Coulomb  collisions  as 


5.(^)  = 


E  L 

"  ^mm  J- 


f{E,pL,s)Ec  OU  idpt , 


(9) 


where  [Timin,  Emax]  is  the  range  of  the  energy  bins  used  in  the 
simulation,  and  the  electron  distribution  function  /(£,  /x,  y)  can 
be  obtained  from  the  corresponding  electron  flux  E{E,  /x,  s)  via 
Equation  (3). 

Alternatively,  one  can  calculate  the  net  downward  energy  flux 
carried  by  the  electrons. 


an 

F.r,is)  = 

a{s) 


P  ^max  P  1 

E  L 


PlEE(E,  fi,  s)dpL , 


(10) 


Se(s)^-dEers{s)/ds,  (11) 

where  fiE F{E ,  pL,  s)  is  the  energy  flux  along  the  loop  and 
the  factor  ao/a(s)  accounts  for  the  variation  of  the  loop 
cross-sectional  area.  This  approach  is  practically  equivalent  to 
Equation  (9),  because  in  the  HXR  energy  range  the  combination 
of  synchrotron  and  bremsstrahlung  radiation  only  constitutes  a 
negligible  fraction  (<10“"^)  of  the  energy  loss  of  a  fast  electron 
due  to  Coulomb  collisions. 

3.2.  Code  Communication 

It  is  desirable  that  the  particle  and  HD  codes  communicate  at 
each  time  step  during  the  time  advance.  The  current  transport 
code,  however,  can  only  provide  a  steady-state  solution  and  does 
not  have  a  time-dependent  capability.  This  can  be  remedied  by 
carefully  selecting  the  communication  time  interval  At,  because 
particle  transport  occurs  on  a  much  shorter  timescale  than 
hydrodynamics.  This  interval  should  be  as  short  as  possible 
provided  that  a  steady-state  transport  solution  can  be  reached. 
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By  considering  the  electron  “lifetime”  (see  Equation  (9)  of 
Petrosian  1973),  which  is  determined  by  the  energy-loss  time 
in  a  given  loop  geometry  and  atmospheric  density  distribution, 
Liu  (2006,  his  Section  7.2.4)  found  the  optimal  interval  to  be 
Af  =  2  s. 

The  remaining  question  is  what  heating  rate  the  HD  code 
should  use  during  its  time  advance  between  adjacent  com¬ 
munications  with  the  particle  code.  Let  us  first  change  the 
independent  variable  of  Se  from  distance  s  to  column  depth 
N(s)  =  /g  ne(s')ds', 

Se[Nis)]^Se{s)/nAs),  (12) 

noting  Seis)ds  =  SeiN)dN  and  dN  =  ngds.  Here  Se{N)  is 
in  units  of  erg  s“'.  We  have  assumed  a  loop  of  uniform  cross- 
section  and  thus  no  magnetic  mirroring,  and  here  we  further 
neglect  synchrotron  loss  Esynch  (valid  for  <1  MeV  electrons). 
Under  these  assumptions,  the  electron  flux  F[E,  /r,  s(N)]  = 
F[E,  fx,  N{s)]  is  a  function  of  column  depth  N,  independent  of 
distance  s,  and  so  does  the  heating  rate  SeiN)  calculated  from 
E{E,ti,N). 

The  HD  response  timescale  is  characterized  by  sound  travel 
time,  which  is  84  s  for  a  sound  speed  of  166  km  s“^  at  T  = 
10^  K  in  a  14  Mm  long  loop  here.  Since  Af  =  2  s  <$C  84  s, 
we  assume  that  Se(N)  is  constant  in  time  between  adjacent  code 
communications.  During  this  Af,  the  spatial  distribution  of  the 
heating  rate  Sds,  t)  varies  with  time  merely  according  to  the 
redistribution  of  density  and  thus  the  variation  of  column  depth, 

Se(s,t)=  Se{N)ne(s,t).  (13) 

In  practice,  at  a  given  time  f  and  distance  s,  we  first  identify  its 
column  depth  N{s,  f),  which  is  then  used  to  look  up  the  heating 
rate  Se(N),  and  then  we  apply  the  local  density  to  obtain  Sds,  t) 
by  Equation  (13). 

Communications  between  the  two  codes  are  summarized  in 
Figure  3.  First,  the  HD  code  passes  the  initial  density  distribution 
to  the  particle  code,  which  then  runs  its  first  steady-state 
calculation  and  returns  the  heating  rate  Se(N)  as  a  function 
of  column  depth  to  the  HD  code.  Next,  the  HD  code  repeatedly 
converts  Se{N)  to  S^is,  t)  as  a  function  of  distance  at  each 
time  step  using  the  latest  density  profile.  Once  the  HD  code 
advances  a  time  interval  of  Af  =  2  s,  it  passes  the  updated 
density  distribution  back  to  the  particle  code,  which  starts  the 
next  cycle  of  iteration. 

4.  SIMULATION  RUNS 

We  have  performed  three  simulation  runs  (see  Table  1)  to 
test  the  relative  effects  of  different  processes.  (1)  In  the  first 
run,  which  we  refer  to  as  the  Old  Model  (abbreviated  by  “O”), 
we  assumed  an  injection  of  electrons  of  a  power  law  at  the 
LT,  and  evaluated  the  heating  rate  and  HD  response  as  in  the 
MEL89  model.  This  model  does  not  calculate  particle  transport 
properly.  (2)  In  the  second  simulation,  we  still  injected  power- 
law  electrons  and  evaluated  electron  transport  and  heating  along 
the  loop  using  our  transport  code.  We  call  this  the  Hybrid 
Model  (or  “H”).  (3)  Finally,  we  employed  our  most  realistic 
model,  where  we  evaluated  from  our  acceleration  code  (PL04) 
the  spectra  of  electrons  at  the  LT  acceleration  site  and  those 
escaping  the  LT  region.  We  also  calculated  transport  and  heating 
using  our  transport  code.  We  call  it  the  New  Model  (or  “N”). 
In  all  cases,  we  assumed  an  identical  initial  HD  state  as  shown 
in  Figure  1 ,  and  calculated  the  HD  response  using  the  MEL89 


Figure  3.  Task  flow  chart  for  communications  between  the  particle  and  HD 
codes. 


code.  We  assumed  the  dynamic  or  modulation  profile  of  the 
number  of  injected  electrons  (power  law  for  Runs  O  and  H  and 
Q  for  Run  N)  to  be  a  triangular  shape  with  a  rise  and  fall  to  be 
30  s  each.  Beyond  this  first  60  s  of  the  impulsive  phase,  while 
the  electron  heating  rate  Se  was  set  to  zero,  we  continued  the 
computation  into  the  decay  phase  until  f  =  90  s. 

4.1.  Run  O:  Old  Model 

We  first  computed  the  HD  response  using  the  same  heat¬ 
ing  function  and  almost  identical  control  parameters  as  the 
“Reference  Calculation”  of  MEL89.  This  model  is  based  on 
the  analytical  solution  of  Emslie  (1978)  that  includes  only  en¬ 
ergy  loss  and  pitch-angle  growth  of  injected  electrons.  By  the 
MEL89  assumption,  the  injected  electron  flux  is  downward 
beamed  (/xo  =  1),  and  its  energy  spectrum  (see  Figure  4)  is 
a  broken  power  law  oc  E~^  with  an  index  5  =  6  above  and 
5  =  —2  below  the  cutoff  (“knee”)  energy  Ee  —  15  keV.  Here 
the  differences  are:  (1)  our  peak  energy  flux,  i.e.,  parameter  U  in 
Equation  (9)  of  MEL89,  is  2.67  x  10'°  erg  cm“^  s“',  while 
they  used  5  x  10^°  erg  cm“^  s“';  and  (2)  we  assume  a  fully 
ionized  hydrogen  plasma  while  they  included  helium  of  6.3%  of 
hydrogen  in  number  density.  The  latter  difference  only  changes 
the  absolute  mass  density  by  1 1%,  while  the  relative  differences 
between  our  models  may  not  be  affected. 

4.2.  Run  H:  Hybrid  Model 

Here  we  used  our  Fokker-Planck  transport  code  in  place  of 
the  approximate  analytical  expression  used  above  to  evaluate 
the  heating  rate.  We  injected  electrons  of  an  identical  power- 
law  spectrum  with  a  narrow-Gaussian  (cr^  =  0.01)  pitch-angle 
distribution  to  emulate  the  beamed  distribution  in  Run  O.  The 
main  difference  here  is  that  the  transport  code  properly  treats 
the  diffusion  process  of  pitch-angle  change  due  to  Coulomb 
collision.  For  the  particle  code,  the  energy  space  is  divided  into 
200  uniform  logarithmic  bins  in  the  range  of  5 1 1  x  10“^-5 1 1  x 
10^^  keV,  while  the  pitch-angle  space  is  divided  into  100  (versus 
24  used  in  Liu  2006)  uniform  bins  in  the  [0,  n]  range. 

4.3.  Run  N:  New  Model 

This  is  a  typical  simulation  using  our  new  model.  It  is  the 
same  as  Run  H,  except  that  the  injected  beamed  power-law 
electron  flux  is  replaced  with  an  isotropic  flux  given  by  the 
stochastic  acceleration  code.  We  used  the  same  acceleration 
parameters  as  PL04  (see  their  Figure  12),  i.e.,  the  characteris¬ 
tic  acceleration  rate  t“'  =  70  s“',  Ug  —  1.5  x  10*°  cm“^, 
Bq  —  400 G,  kt,T  —  1.53  keV,  and  the  acceleration  re¬ 
gion  size  L  =  5  X  10**  cm.  We  modulated  the  rate  (Q,  see 
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Table  1 

Summary  of  Simulation  Runs 


Runs 

Injected  Electron 

Particle 

^max 

^min 

fD>100 

^apex 

2max 

^e.max 

Spectrum 

Ang.  Distr. 

Transport 

(km  s-l) 

(S) 

(km  s-l) 

(s) 

(S) 

(10°  K) 

(10'“  cm-3) 

0(ld) 

Power  law 

Beamed 

Analy.  Approx. 

565 

35 

-115 

to 

29 

21.1 

6.96 

H(ybrid) 

Power  law 

Beamed 

Fokker-Planck 

570 

35 

-106 

9 

28 

22.0 

7.44 

N(ew) 

Stock.  Accel. 

Isotropic 

Fokker-Planck 

598 

36 

-102 

9 

23 

25.9 

8.44 

Notes.  Runs  O  and  H  have  an  identical  injected  power-law  electron  flux,  with  a  spectral  index  8  =  6  and  low-energy  cutoff  Ec  =  15  keV;  Umax 
and  tvraax'  iTiaximum  (upflow,  u  >  0)  velocity  and  its  time  stamp;  Umm^  minimum  (downflow,  u  <  0)  velocity  (in  the  upper  chromosphere);  /y>ioo- 
time  when  the  upflow  velocity  exceeds  100km  s“^;  tapex-  time  when  the  evaporation  front  (density  jump)  reaches  the  loop  apex;  Tmax  and  max- 
maximum  coronal  temperature  and  electron  density.  All  runs  have  the  same  peak  energy  deposition  (electron  heating)  flux  for  the  loop  as  a  whole, 
f^max  =  2.67  X  10^®  erg  s“^  cm“^. 
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Figure  4.  Electron  flux  spectra  F{E)  times  i?'  at  the  peak  time  (r  =  30  s): 
(1)  acceleration  region  flux  Fac.  escaping  flux  Fesc,  and  equivalent  thick-target 
flux  Fthick  for  Run  N,  and  (2)  injected  flux  for  Runs  O  and  H  (5  =  6  above 
Ec  =  15  keV).  The  triple-dot-dashed  line  (arbitrary  scale)  represents  the  source 
term  Q(E)  of  a  Maxwellian  (thermal,  kjF  =  1.53  keV)  distribution  used  in 
Run  N,  which  peaks  at  F  =  3f:i,T  in  this  E^E{E)  plot. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Equation  (1))  of  electrons  being  supplied  to  the  acceleration 
region  with  the  same  triangular  time  profile,  such  that  the  peak 
electron  heating  flux  Emax  equals  that  of  Run  O  or  H. 

Figure  4  shows  various  electron  flux  spectra  used  in  this  study. 
In  comparison  with  the  background  thermal  distribution,  both 
the  acceleration  region  flux  Fac  and  escaping  flux  Fesc  have  a 
quasi-thermal  component  that  smoothly  extends  to  a  nonthermal 
tail  at  high  energies.  Fesc  is  smaller  than  Fac  and  their  relative 
difference  decreases  with  energy  due  to  the  energy-dependent 
confinement  of  electrons  by  turbulence  in  the  acceleration  region 
(see  Equation  (2)).  Unlike  that  (dot-dashed)  in  Run  O  or  H, 
the  flux  Fesc  injected  into  the  transport  region  does  not  invoke 
any  arbitrary  low-energy  cutoff.  The  two  fluxes,  however,  have 
similar  slopes  in  the  intermediate  energy  range  around  20  keV. 
The  equivalent  thick-target  electron  flux  Fthkk  (see  Section  2.4), 
as  expected,  has  a  harder  spectrum  than  Fac  and  Fesc  in  the 
10-1000  keV  range. 

5.  SIMULATION  RESULTS:  COMPARISON  OF  FLUID 
DYNAMICS 

To  determine  how  much  our  proper  transport  and  acceleration 
calculations  affect  the  atmospheric  response,  we  compare  the 


s  (Mm)  s  (Mm)  s  (Mm)  s  (Mm) 


Figure  5.  Evolution  of  plasma  electron  density,  temperature,  pressure,  velocity, 
and  electron  heating  rate  for  Runs  H  and  N. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


results  of  the  three  simulations,  using  Run  H  as  the  reference 
case. 


5.1.  Hydrodynamic  Evolution 

As  shown  in  Figure  5,  Run  H  exhibits  similar  general  HD 
evolution  as  described  in  MEL89  (their  Figures  1  and  2). 
Electron  heating  {Se)  is  initially  concentrated  in  the  upper 
chromosphere,  producing  overpressure  that  drives  an  upflow 
(u  >  0)  and  a  recoil  downflow  (v  <  0).  At  f  =  9  s,  the 
upflow  velocity  exceeds  100  km  s“'  and  a  density  jump  or 
evaporation  front  has  developed  slightly  above  the  TR.  It  travels 
upward  and  reaches  the  loop  apex  at  f  =  28  s.  The  density 
jump  is  then  reflected  back  and  the  material  piles  up  due 
to  the  reflective  boundary  condition  imposed,  which  can  be 
understood  as  plasma  flow  from  the  other  half  of  the  loop 
in  the  assumed  symmetric  geometry.  The  upflow  reaches  its 
maximum  velocity  of  u^ax  =  570  km  s“i  at  r  =  35  s,  which 
is  delayed  by  5  s  from  the  energy  deposition  peak  at  f  =  30  s. 
Chromospheric  evaporation  then  gradually  subsides.  These 
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Figure  6.  History  of  plasma  electron  density,  temperature  and  velocity  at  5  =  1 
Mm  from  the  injection  site  for  different  simulations. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


features  of  the  temporal  evolution  can  also  be  seen  from 
the  history  of  various  quantities  at  a  hxed  position  in  the 
upper  corona  as  shown  in  Figure  6.  Note  that,  late  during  the 
simulation,  the  coronal  temperature  gradually  decreases  mainly 
due  to  conductive  cooling,  while  the  coronal  density  continues 
to  increase,  even  after  the  cease  of  electron  heating  at  f  =  60  s. 
This  is  caused  by  sustained  chromospheric  evaporation  that 
results  from  heating  of  the  chromosphere  by  the  same  conductive 
flux  that  cools  the  hot  corona. 

In  comparison,  we  find  that  Run  O  is  very  similar  but  less 
intense  than  to  Run  H  (see  Figure  6  and  Table  1).  This  can  be 
more  clearly  seen  from  various  HD  quantities  and  heating  and 
cooling  rates  during  the  early  phase  shown  in  Figure  7.  The 
small  differences  (on  the  order  of  10%)  between  the  two  cases 
are  mainly  caused  by  the  slight  overestimate  of  electron  heating 
of  Run  O  in  the  chromosphere  (Figure  7(d)),  due  to  its  inaccurate 
way  of  calculating  particle  transport  noted  earlier.  This  indicates 
that,  for  this  specific  case,  the  analytical  heating  function  used 
in  MEL89  provides  an  acceptable  approximation  for  the  more 
accurate  Fokker-Planck  calculation. 

In  contrast,  the  HD  evolution  of  Run  N  is  faster  and  more 
intense  than  that  of  Run  H  (Figures  5  and  6,  Table  1),  despite 
the  same  peak  value  of  the  spatially  integrated  electron  heating 
rate  for  the  two  cases.  These  differences  are  primarily  due  to  the 
different  energy  spectra  of  injected  electrons,  further  enhanced 
by  their  different  pitch-angle  distributions. 

In  Run  N,  electrons  of  a  few  keV  in  the  quasi-thermal 
component  of  the  spectrum  (Figure  4)  are  small  in  energy  but 
large  in  number  and  thus  dominate  the  total  energy  content. 
These  electrons  produce  heating  at  relatively  small  column 
depths  by  Coulomb  collisions.  This  occurs  high  in  the  corona 
(Figure  7(d)),  where  the  radiative  loss  rate  Trad  (Figure  7(e))  is 
relatively  small  due  to  the  low  density  and  high  temperature. 
As  a  result,  significant  net  heating  sets  in  there,  which  leads 
to  a  local  temperature  (Figure  7(b))  and  pressure  surge.  This 


local  coronal  heating  is  enhanced  by  the  large  effective  column 
depth  (Neff  —  N /(pi),  where  (/r)  is  the  mean  pitch-angle 
cosine)  resulting  from  the  isotropic  pitch-angle  distribution  of 
the  injected  electrons.  The  increased  temperature  leads  to  a  large 
downward  heat  conduction  flux  (Figure  7(f)),  and  the  pressure 
gradient  force  drives  a  downward  mass  flow  in  the  high  corona 
(Figure  7(c)). 

In  Run  H,  contrastingly,  the  electron  spectrum  (Figure  4) 
peaks  at  the  cutoff  energy  Ec,  which  leads  to  a  deficit  in  low- 
energy  electrons.  In  addition,  the  pitch-angle  distribution  here 
is  beamed  (rather  than  isotropic).  This  electron  population,  on 
average,  penetrates  deeper  into  the  atmosphere  than  that  in  Run 
N  and  deposits  its  energy  primarily  in  the  upper  chromosphere. 
This  results  in  less  heating  in  the  corona  and  stronger  and  more 
widespread  heating  in  the  chromosphere  (Figures  7(b)  and  7(d)). 
Consequently,  in  spite  of  the  larger  and  broader  radiative  cooling 
(Figure  7(e)),  the  local  overpressure  in  the  chromosphere  is 
stronger  than  that  in  Run  N  early  on,  which  drives  a  higher 
velocity  upflow  (Figure  7(c),  r  =  6  s).  Also,  unlike  Run  N,  there 
is  no  significant  downward  coronal  heat  conduction  (Figure  7(f)) 
or  mass  flow. 

5.2.  Heating  and  Cooling 

A  remaining  question  is  why  Run  N  has  more  dramatic  overall 
HD  changes  in  the  long  run.  To  answer  this,  we  examine  the 
relationship  between  different  energy  gain  and  loss  terms — 
electron  heating,  radiative  loss,  and  conductive  heating  and 
cooling — particularly  early  in  the  flare  and  near  the  TR  where 
chromospheric  evaporation  takes  place. 

InRunNatr=  1  s  (see  top  panel  ofFigure  8(a)),  the  electron 
heating  rate  Se  peaks  in  the  TR  because  of  the  sharp  increase 
there  in  ambient  density  and  associated  collisional  energy  loss  of 
energetic  electrons.  So  does  the  radiative  loss  rate  Trad,  since  it  is 
proportional  to  n^np  and  0(7’)  that  peaks  at  T  ~  (1-3)  x  10^  K 
which  is  the  TR  temperature.  However,  due  to  their  different 
functional  dependencies  on  density  and  temperature,  Se  peaks 
at  a  slightly  lower  position  than  Their  combination  Se  —  Trad 
(panel  2,  Figure  8(a))  thus  results  in  cooling  in  the  upper  TR  and 
heating  in  a  shallow  layer  below  it  in  the  upper  chromosphere. 
Meanwhile,  the  conductive  flux  carries  energy  from  the  hot 
upper  corona  to  the  upper  TR  where  localized  heating  (Scond) 
is  produced  and  counteracts  radiative  cooling.  (Conduction 
is  prohibited  in  the  chromosphere  where  the  temperature  is 
maintained  at  10^  K.)  The  net  energy  gain  resulting  from 
the  interplay  of  electron  heating,  radiative  cooling,  and  heat 
conduction  is  thus  localized  in  the  upper  chromosphere,  where 
temperature  is  raised  substantially  (panel  3,  Figure  8(a)).  This 
leads  to  a  local  pressure  hump  which  drives  an  upflow  into 
the  corona  and  a  downflow  into  the  chromosphere  (panel  4, 
Figure  8(a)). 

As  time  proceeds  (see  Figure  8(b))  and  chromospheric  ma¬ 
terial  is  being  heated  from  T  =  10^  K  to  ~10^  K  where  0(7’) 
reaches  its  maximum,  radiative  loss  gradually  overtakes  elec¬ 
tron  heating  in  the  TR  and  upper  chromosphere.  This  means  that 
energy  directly  deposited  by  electrons  in  these  places  is  imme¬ 
diately  radiated  away  and  very  little  is  left  to  heat  the  plasma. 
However,  as  we  noted  above,  a  significant  portion  of  the  energy 
content  of  the  injected  electrons  in  Run  N  is  deposited  in  the 
upper  corona  (Figure  7(d))  where  radiative  loss  is  negligible  and 
then  transported  by  heat  conduction  to  the  lower  atmosphere.  In 
time,  conduction  plays  an  increasingly  important  role  in  heat¬ 
ing  the  lower  corona  and  TR  as  it  gradually  exceeds  the  net 
direct  heating  or  combined  electron  heating  and  radiative  loss 
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Figure  7.  Same  as  Figure  6  but  for  snapshots  of  various  quantities  as  a  function  of  distance  at  selected  times  early  in  the  flare. 
(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Se  —  Erad  (Eigure  8(b))  in  these  regions.  Because  of  this,  as  of 
t  —  1  s  the  location  of  the  primary  net  heating,  i.e.,  the  highest 
local  overpressure,  and  the  maximum  downflow  velocity  have 
shifted  from  the  upper  chromosphere  up  to  the  TR.  In  fact,  the 
peak  of  the  conductive  flux  Fcond  (see  Figure  7(f))  and  the  re¬ 
gion  of  significant  conductive  heating  (Scond  oc  —dF^omi/ds) 
below  it  are  initially  located  in  the  upper  corona  and  propagate 
down  the  loop  with  time.  As  the  Fcond  peak  approaches  the  TR 
at  ~10  Mm  during  the  interval  of  8-10  s,  the  maximum  up- 
flow  velocity  Umax  in  the  loop  increases  abruptly  (see  Figure  9). 
This  further  emphasizes  the  role  played  by  heat  conduction  here 
in  redistributing  energy  deposited  by  electrons  and  in  driving 
chromospheric  evaporation. 

In  Run  H,  the  injected  electron  flux  with  a  low-energy  cutoff 
has  profound  consequences.  As  noted  earlier,  lack  of  low-energy 
electrons  makes  the  chromosphere,  rather  than  the  upper  corona, 
the  primary  location  of  direct  heating  early  in  the  flare  (see 
Figure  7(d)).  Net  direct  heating  Se  —  and  the  increase  of 
local  temperature  and  pressure  extend  from  the  TR  to  much 
deeper  layers  of  the  chromosphere  than  in  Run  N  (Figure  8(c)). 
Moreover,  since  the  coronal  temperature  does  not  increase 
rapidly,  the  downward  conductive  flux  here  is  more  than  an  order 
of  magnitude  smaller  (Figure  7(f)).  Net  direct  heating  generally 
dominates  over  conductive  heating  when  integrating  over  the 


volume  of  the  lower  atmosphere.  As  a  result,  in  comparison 
with  Run  N,  a  relatively  larger  portion  of  the  total  energy 
content  of  the  injected  electrons  is  lost  in  radiative  cooling. 
This  is  why  the  overall  HD  development  here  is  more  gentle 
than  that  of  Run  N,  despite  the  fact  that  they  have  the  same 
energy  deposition  flux.  The  primary  underlying  physics  is  their 
different  spatial  distributions  of  electron  heating  Se  caused  by 
their  different  electron  injections.  Note  that  MEL89  obtained 
qualitatively  similar  results  when  different  values  of  the  cutoff 
energy  or  spectral  index  were  considered. 

We  note  that,  later  during  the  impulsive  phase  in  Run  H 
(Figure  5),  as  the  coronal  density  has  increased  considerably  due 
to  chromospheric  evaporation,  relatively  more  electron  energy 
is  directly  dumped  in  the  corona  while  less  in  the  chromosphere. 
Consequently,  conductive  heating  in  the  TR  becomes  important. 
In  this  sense,  the  physical  distinction  between  the  two  runs 
gradually  diminishes  in  the  late  stage  of  the  flare. 

5.3.  Velocity  Distribution 

Here  we  examine  observables  that  can  be  checked  against 
data:  (1)  the  temperature  dependence  of  the  plasma  velocity 
(Figure  10,  left);  and  (2)  the  velocity  differential  emission 
measure  (VDEM;  Figure  10,  right)  defined  by  Newton  et  al. 
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(a)  Run  N  (t=1  s)  (c)  Run  H  (t=l  s) 


s  (Mm)  s  (Mm) 

Figure  8.  Comparison  between  Runs  N  (left)  and  H  (right)  of  a  detailed  view 
near  the  TR  at  two  selected  times.  At  each  time,  the  quantities  in  the  four  panels 
(from  top  to  bottom)  are:  (1)  electron  heating  rate  Se  and  radiative  loss  rate  Lrad; 
(2)  Se  —  Lrad  (dotted)  and  conductive  heating  rate  5cond  (solid);  (3)  electron 
density  (left  scale)  and  temperature  (right  scale);  and  (4)  pressure  (left)  and 
velocity  (right).  Note  different  vertical  scales. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 

(1995)  as  the  emissivity  (G[7’])  weighted  measure  of  material 
with  line-of-sight  (LOS)  velocity  Vlos-  Assuming  that  the 
flare  is  located  at  the  solar  disk  center  such  that  the  LOS 
is  perpendicular  to  the  loop  apex,  we  calculated  the  specific 
VDEM/a(i')  for  the  Ca  xix  3.18  A  line  following  Newton  et  al. 
Here  G{T)  is  given  by  the  Chianti  package  (Young  et  al.  2003). 

Early  in  the  flare  (Eigures  10(a)  and  10(c)),  Run  H  has 
higher  upflow  velocities  and  downflow  temperatures  than  Run 
N,  owing  to  its  stronger  net  direct  heating  with  a  deeper  extent 
in  the  chromosphere  (Figure  8(c)).  Run  N  has  an  additional  hot 
downflow  component  at  >2  MK,  due  to  the  expansion  of  the 
heated  upper  corona  mentioned  earlier.  This  hot  component  also 
dominates  the  VDEM  at  Ulos  <  0  (Figure  10(c))  because  of  the 


Time  (s) 

Figure  9.  Early  history  of  the  maximum  upflow  velocity  Umax  (left  scale)  and 
the  position  of  the  maximum  conductive  flux  (right  scale)  for  Run  N.  The  two 
vertical  lines  mark  the  inteiwal  when  Umax  experiences  a  sharp  increase  while 
the  maximum  Fcond  rapidly  approaches  the  TR  at  5  =  ~  10  Mm. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 

sharp  rise  of  G{T)  up  to  T  =  29  MK.  Later  at  t  =  28  s 
(Figure  10(b)),  Run  N  overtakes  H  in  upflow  velocity  and 
temperate,  while  its  hot  coronal  downflow  has  disappeared.  Run 
N  has  a  bimodal  VDEM  (Figure  10(d))  with  a  strong  stationary, 
hot  component  located  near  the  loop  apex  (also  see  Figure  5), 
while  Run  H  exhibits  a  more  gradual  progression  toward  high 
velocities.  This  distinction  is  similar  to  what  Newton  et  al.  (1995, 
see  their  Figure  2)  found  for  models  with  different  initial  coronal 
densities. 

6.  EFFECTS  OF  HYDRODYNAMICS  ON  PARTICLE 

TRANSPORT  AND  X-RAY  EMISSION 

We  now  turn  our  attention  to  the  effects  of  fluid  dynamics 
on  particle  characteristics,  namely,  electron  transport  and  non- 
thermal  radiation.  Here  we  present  the  result  of  Run  N  as  a 
typical  example.  We  will  examine  first  the  energy  distributions 
of  electrons  and  bremsstrahlung  photons,  and  then  their  spatial 
distributions. 

6.1.  Electron  and  Photon  Energy  Spectra 

Figure  11  (top  panels)  shows  the  evolution  of  the  angle- 
integrated  electron  flux  spectrum  E^E(E,  s)  at  different 
locations  in  the  loop.  In  general,  at  a  given  time  and  a  mod¬ 
erately  large  distance  or  column  depth,  there  is  a  deficit  in  low- 
energy  electrons  due  to  collisional  energy  losses  and  scatterings 
on  the  paths  of  the  electrons  from  the  injection  site.  This  ap¬ 
pears  as  a  turnover  in  the  spectrum  and  slight  spectral  hardening 
just  above  the  turnover  energy.  As  distance  increases,  progres¬ 
sively  more  low-energy  electrons  are  lost,  and  thus  the  overall 
flux  decreases  while  the  spectral  turnover  shifts  to  higher  en¬ 
ergies  (Leach  &  Petrosian  1981).  At  f  =  2  s  (Figure  11(a)) 
the  TR  is  located  at  distance  Str  =  9.97  Mm  and  column  depth 
Wr  =  6.1  X  10'*  cm“^.  The  two  fluxes  di  s  —  A  and  10  Mm 
(thin  solid,  dotted)  are  located  in  the  low-density  corona  or  TR 
at  small  column  depths  from  the  injection  site,  and  thus  appear 
similar  in  shape  to  the  escaping  flux  Tgsc  (thick  dotted).  Other 
fluxes  at  s  =  11, 12,  and  1 3  Mm  are  located  in  the  chromosphere 
at  large  column  depths,  and  thus  exhibit  substantial  reduction  of 
low-energy  electrons. 

In  this  study,  the  flux  (Uesc)  injected  into  the  transport  region 
does  not  change  with  time  in  spectral  shape  but  only  varies 
in  normalization.  So  does  the  electron  flux  at  a  given  column 
depth.  However,  as  chromospheric  evaporation  develops,  the  TR 
retreats  to  lower  altitudes  (Figure  5),  while  the  coronal  density 
increases.  This  causes  the  change  of  the  column  depth  and  thus 
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Figure  10.  Comparison  between  Runs  H  and  N:  velocity  vs.  temperature  (left)  and  corresponding  specific  velocity  differential  emission  measure  (photons 
cm~^  s~^  sr“^  [km  vs.  LOS  velocity  (right). 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Figure  11.  Spectra  of  angle-integrated  electron  flux  (top)  and  con'esponding  photon  emission  rate  (bottom)  multiplied  by  energy  squared  at  three  selected  times 
for  Run  N.  The  thick  dashed  line  represents  the  LT  (acceleration  region)  spectrum  (Fac  or  /lt),  while  the  thin  colored  lines  (solid,  dotted,  dashed,  dot-dashed,  and 
triple-dot-dashed)  are  the  spectra  at  distances  s  =4,  10,  11,  12,  and  13  Mm  from  the  injection  site.  The  thick  dotted  line  indicates  the  escaping  electron  flux  (Fesc)  in 
the  top  panels  and  the  equivalent  FP  photon  emission  rate  (/pp)  in  the  bottom  panels.  The  legends  show  the  current  values  of  the  position  (%)  and  column  depth  (Atr) 
of  the  TR. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


the  electron  spectrum  at  a  fixed  location,  which  is  what  we  notice 
here.  At  t  =  30  s  (Eigure  11(b))  when  the  TR  shifts  down  to 
Jtr  =11.3  Mm  at  =  4.2  x  10'®  cm“^,  the  spectrum  at  = 
1 1  Mm  looks  similar  to  the  other  coronal  spectra  (at  i  =  4 
and  10  Mm)  but  different  from  the  chromospheric  spectra. 
Meanwhile,  the  relative  difference  between  the  first  coronal 


spectrum  (s  =  4  Mm)  and  the  injected  spectrum  F^sc  becomes 
larger,  because  of  the  increasing  coronal  density  and  column 
depth  between  them.  This  trend  continues  through  the  end  of 
the  simulation. 

The  bottom  panels  of  Eigure  11  show  the  correspond¬ 
ing  bremsstrahlung  photon  spectra  e^I(e,  s)  defined  in 
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Figure  12.  Spatial  distributions  of  angle-integrated  electron  flux  (top)  and  photon  emission  rate  (bottom)  at  thi'ee  selected  times  for  Run  N.  From  top  to  bottom, 
different  lines  in  each  panel  represent  electron  or  photon  energies  of  3.1,  6.1,  12.3,  24.5,  48.8,  97.4,  and  294.1  keV.  The  step  in  the  5  <  0  region  corresponds  to  one 
half  length  of  the  acceleration  region  at  the  top  of  the  loop  (see  Figure  1).  The  step  at  the  TR  is  due  to  the  jump  in  the  ambient  density. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Equation  (7).  Like  its  electron  counterpart,  the  LT  photon  spec¬ 
trum  €^/lt  (thick  dashed)  shows  a  thermal-like  component  at 
low  energies  and  a  nonthermal  tail  at  intermediate  energies  with 
a  cutoff  at  high  energies.  The  equivalent  FP  spectrum  e^/pp  be¬ 
low  Sti  (thick  dotted,  defined  in  Section  2.4)  is  harder  than  the 
LT  spectrum.  Their  shapes  in  the  tens  to  hundreds  keV  range  are 
commonly  observed  for  LT  and  FP  sources.  At  larger  distances, 
the  spectra  (thin  lines)  become  progressively  harder  because  the 
corresponding  electron  spectra  have  the  same  variation  (Leach 
&  Petrosian  1983)  as  we  have  noticed  above,  and  high-energy 
(>  tens  keV)  emission  increases  because  of  higher  ambient  den¬ 
sities.  As  the  coronal  density  increases  with  time,  for  the  same 
reason  mentioned  above,  photon  spectra  at  positions  originally 
located  in  the  corona  become  progressively  harder.  At  the  same 
time,  as  the  TR  treats,  spectra  at  positions  originally  in  the  chro¬ 
mosphere  but  now  in  the  corona  become  softer.  In  other  words, 
the  difference  in  spectral  shape  between  X-ray  emissions  at  dif¬ 
ferent  positions  diminishes  with  time  (Figures  1 1(e)  and  1 1(f)). 

6.2.  Electron  and  Photon  Spatial  Distributions 

Figure  12  (top  panels)  shows  the  evolution  of  the  angle- 
integrated  electron  flux  F(E,  s)  as  a  function  of  distance  at  dif¬ 
ferent  energies  from  3  to  300  keV.  The  step  at  5  =  0  is  owing  to 
the  difference  between  the  LT  (acceleration  region)  flux  Fac  and 
the  escaping  flux  Fgsc  mentioned  before  (see  Figure  4).  In  gen¬ 
eral,  the  electron  flux  decreases  with  distance  or  column  depth 
from  the  injection  site.  The  slope  of  the  curve  —dF{E,  s)/ds 
is  steeper  at  lower  energies  because  lower-energy  electrons  lose 
energy  faster  and  are  more  sensitive  to  pitch-angle  scattering. 
At  a  given  energy,  the  slope  depends  on  the  ambient  density, 
because  dF(E,  s)/ds  —  ne{dF{E,  N)/dN'\,  where  F{E,  N)  is 
generally  a  smooth  function  of  N  (McTieman  &  Petrosian  1 990). 
Therefore,  if  there  is  a  rapid  increase  in  density  with  distance, 
the  slope  increases  quickly;  the  opposite  would  be  true  if  the 


density  were  to  decrease  sharply.  A  sharp  break  is  obvious  at 
the  TR,  and  a  milder  break  occurs  at  the  evaporation  front  (e.g., 
at  s  ~  5  Mm  in  Figure  12(b)).  An  upward  break  or  flattening 
is  visible  at  ~  3  Mm  where  the  density  jump  is  reflected 
back  from  the  loop  apex  at  t  =  30  s  (Figure  12(c)).  It  is  also 
visible  just  below  the  TR  (see  Figures  12(b)  and  12(c)),  where 
a  sharp  density  decrease  from  the  narrow  density  spike  occurs 
(Figure  5).  As  density  increases  due  to  chromospheric  evapo¬ 
ration,  the  spatial  distribution  in  the  corona  becomes  steeper. 
The  slope  variation  with  distance  (i.e.,  breaks)  and  time  is  more 
pronounced  at  lower  energies  for  the  same  reason  noted  above. 

Figure  12  (bottom  panels)  shows  the  evolution  of  the  corre¬ 
sponding  spatial  distribution  of  angle-integrated  bremsstrahlung 
emission  I{e,  s).  In  the  early  stage,  low-energy  emission  comes 
primarily  from  the  LT,  while  high-energy  emission  is  concen¬ 
trated  below  the  TR.  Because  nonthermal  bremsstrahlung  radi¬ 
ation  is  proportional  to  both  the  electron  flux  and  local  target 
density,  photon  emission  reflects  details  of  the  density  distribu¬ 
tion  in  a  more  pronounced  way  than  the  electron  flux  profile. 
As  is  evident,  the  emission  profile  clearly  indicates  density  fea¬ 
tures,  including  the  evaporation  front  and  the  density  spike  just 
below  the  TR.  The  local  emission  enhancement  at  the  evapora¬ 
tion  front  moves  upward  with  time,  which  may  be  responsible 
for  the  observed  X-ray  sources  moving  along  the  loop  (Liu  et  al. 
2006a;  Sui  et  al.  2006)  mentioned  in  Section  1.1.  As  time  pro¬ 
ceeds,  relatively  more  emission  comes  from  the  coronal  portion 
of  the  loop  because  of  the  increased  density.  Specifically,  at 
low  energies,  the  emission  intensity  decreases  with  distance  in 
the  corona  more  sharply  than  before.  At  intermediate  energies, 
we  find  a  femporal  fransition  from  FP-dominated  emission  to 
LT-dominated  emission,  which  occurs  at  progressively  higher 
energies,  reminiscent  of  the  phenomenon  observed  by  Liu  et  al. 
(2006a).  At  high  energies,  such  a  change  is  not  visible  because 
the  high-density  corona  still  looks  transparent  to  high-energy 
electrons,  but  the  retreat  of  the  TR  is  obvious. 
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7.  SUMMARY  AND  DISCUSSION 

We  have  performed  simulations  of  solar  flares  that  self- 
consistently  combine  acceleration,  transport,  and  radiation  of 
energetic  electrons  (using  the  Stanford  unified  code)  with  fluid 
dynamics  of  the  atmospheric  response  (using  the  NRL  flux  tube 
code).  As  the  first  successful  one  of  its  kind,  this  model  improves 
on  previous  HD  simulations  in  two  major  aspects.  First,  it 
includes  more  accurate  evaluation  of  electron  heating  from 
full  Fokker-Planck  calculation  of  particle  transport.  Second, 
it  uses  a  more  realistic  electron  spectrum  from  the  stochastic 
acceleration  model  (PL04)  as  the  injection  to  the  transport 
region.  We  compare  this  more  advanced  treatment  with  models 
in  which  an  ad  hoc  electron  distribution  of  a  power  law  with  a 
low-energy  cutoff  is  injected  into  the  loop  and/or  transport  is 
dealt  with  approximately.  Our  conclusions  are  as  follows. 

1 .  For  the  specific  injection  of  beamed,  power-law  electrons, 
the  old  analytical  model  of  MEL89  (Run  O)  provides  an 
acceptable  approximation.  Its  result  differs  by  ~I0%  from 
that  of  the  reference  hybrid  model  (Run  H)  obtained  by  the 
more  accurate  Fokker-Planck  calculation  (see  Table  1  and 
Figures  6  and  7). 

2.  In  the  new  model  (Run  N),  where  the  injected  electron 
spectrum  is  based  on  stochastic  acceleration,  we  find  higher 
coronal  temperatures  and  densities,  larger  upflow  velocities, 
and  faster  increases  of  these  quantities  than  the  hybrid 
model  (Run  H,  Figure  6).  This  is  mainly  because  the  new 
injected  electron  spectrum  smoothly  spans  from  a  quasi- 
thermal  component  to  a  nonthermal  tail  (Figure  4).  The 
low-energy  electrons  in  the  quasi-thermal  regime,  which 
contain  the  bulk  of  the  total  energy  budget,  deposit  their 
energy  primarily  in  the  corona.  This  results  in  significant 
coronal  heating  and  thus  a  large  downward  heat-conduction 
flux  that  helps  drive  “evaporation”  of  plasma  at  the  TR.  In 
contrast,  the  electron  spectrum  in  the  hybrid  model  with  a 
low-energy  cutoff  leads  to  more  energy  directly  deposited  in 
the  chromosphere,  which  is  quickly  radiated  away,  leaving 
less  energy  to  produce  actual  heating  (Figures  7  and  8).  This 
is  qualitatively  consistent  with  the  conclusion  of  MEL89 
where  an  electron  spectrum  with  a  smaller  low-energy 
cutoff  or  steeper  slope  resulted  in  a  stronger  chromospheric 
evaporation. 

3.  The  energy  and  spatial  distributions  of  energetic  electrons 
and  bremsstrahlung  photons  bear  the  fingerprint  of  the 
changing  density  distribution  caused  by  chromospheric 
evaporation.  In  general,  as  time  proceeds,  the  electron 
and  photon  spectra  at  positions  remaining  in  the  corona 
become  progressively  harder  because  of  the  increasing 
coronal  density,  while  those  at  positions  previously  in  the 
chromosphere  and  now  in  the  corona  (due  to  the  retreat  of 
the  TR)  become  softer  (Figure  11).  Any  density  jump  in 
space  results  in  a  sudden  change  in  the  spatial  distributions 
of  energetic  electrons  and  X-ray  photons  (Figure  12).  In 
particular,  the  evaporation  front  appears  as  a  local  emission 
enhancement,  which,  in  principle,  can  be  imaged  by  X-ray 
telescopes. 

7.7.  Comparison  with  Observations 

Over  several  decades  (see  review  by  Antonucci  et  al.  1999), 
Doppler  observations  have  indicated  hot,  fast  (>100  km  s“') 
upflows  (Doschek  et  al.  1980)  and  cool,  slow  (>10  km  s“') 
downflows  (Wuelser  et  al.  1994;  Brosius  &  Phillips  2004)  dur¬ 
ing  flares.  This  is  consistent  with  chromospheric  evaporation 


and  momentum  recoil  shown  in  our  and  earlier  simulations  (e.g., 
Mariska  et  al.  1982;  Fisher  et  al.  1985c).  In  addition,  Milligan 
&  Dennis  (2009)  reported  plasma  velocities  at  multiple  temper¬ 
atures  obtained  from  Hinode  EUV  Imaging  Spectrometer  (FIS) 
observations,  which  show  excellent  agreement  with  the  Run  N 
curve  in  Figure  10(b)  throughout  the  10^-10^  K  range.  Such  an 
agreement  with  HD  simulations  has  not  been  seen  before. 

Heat  conduction,  compared  with  more  popular  direct  electron 
heating,  plays  an  important  role  in  driving  chromospheric 
evaporation  in  our  new  model.  Observational  support  of  this  was 
first  reported  by  Zarro  &  Lemen  (1988),  who  found  an  upflow 
velocity  <50  km  s“'  at  T  ~  6  MK  during  the  cooling  phase 
of  a  flare.  Milligan  (2008)  recently  reported  an  unusually  high 
temperature  (2  MK)  downflow  at  ~14  km  s“'  in  a  microflare 
with  no  detection  of  HXRs  (implying  a  low  flux  of  nonthermal 
electrons),  which  possibly  results  from  conductive  heating.  This 
downflow  could  be  due  to  the  thermal  expansion  early  in  the 
corona  (Figure  10(a))  or  the  momentum  recoil  later  in  the 
chromosphere  (Figure  10(b)).  More  recently,  Battaglia  et  al. 
(2009)  interpreted  the  growths  of  the  SXR  emission  measure 
observed  early  during  flares  (before  HXRs  being  observed)  as 
new  evidence  of  conduction-driven  chromospheric  evaporation. 

Our  simulations  predict  that  evaporation  upflows  tend  to  have 
higher  temperatures  when  conductive  heating  dominates  over 
direct  electron  heating,  while  the  opposite  is  true  for  recoil 
downflows  (Figure  10).  This  can  be  checked  against  observed 
distributions  of  the  plasma  bulk  velocity  versus  temperature 
(e.g.,  Milligan  &  Dennis  2009).  Our  simulated  VDEM  can 
also  be  readily  used  to  synthesize  emission  lines  (Newton  et  al. 
1995)  to  be  compared  with  observations  and  help  differentiate 
theoretical  models,  because  of  the  sensitive  dependence  of 
VDEM  on  heating  mechanisms. 

7.2.  Future  Work 

This  paper  is  the  first  in  a  series  and  we  have  established 
the  numerical  model  and  presented  initial  results.  In  followup 
papers,  we  will  explore  the  parameter  space  and  use  this  model 
to  investigate  the  Neupert  effect  and  the  observed  moving  X-ray 
sources  (Liu  et  al.  2006a;  Sui  et  al.  2006).  More  importantly,  we 
will  incorporate  the  atmospheric  feedback  on  the  acceleration 
process.  This  is  because  chromospheric  evaporation  may  change 
the  physical  condition  (e.g.,  plasma  density  and  temperature) 
in  the  LT  acceleration  region.  The  density  enhancement,  for 
example,  causes  the  ratio  of  electron  plasma  frequency  to  gyro- 
frequency  a  =  cOpg/Q.^  oc  /Bq  to  increase.  This  can  lead  to 
the  reduction  of  the  efficiency  of  electron  acceleration  (PL04) 
and  thus  the  quenching  or  spectral  softening  (e.g.,  Parks  & 
Winckler  1969)  of  nonthermal  HXR  tails  observed  during  the 
late  stages  of  flares. 

Some  technical  aspects  of  this  model  can  be  improved  in  the 
future.  (1)  The  fully  ionized  hydrogen  plasma  assumed  here 
can  be  replaced  by  a  plasma  of  a  more  realistic  solar  abundance, 
with  the  inclusion  of  neutrals  and  ionization  equilibrium.  (2)  The 
“cold”  target  assumption  in  the  transport  code  can  be  abandoned 
and  replaced  with  Coulomb  diffusion  in  energy  space  (Spitzer 
1962;  Miller  et  al.  1996)  for  a  general  “warm”  target  plasma 
(Emslie  2003).  (3)  In  particle  transport'^  and  HD  calculations, 
one  can  include  energetic  protons  (compare,  Emslie  et  al.  1998) 
and  heavier  ions,  whose  momentum  loss  to  the  background 
plasma,  in  addition  to  the  overpressure  (Kosovichev  2006) 

^  Stochastic  acceleration  of  ions  has  been  modeled  (Miller  &  Roberts  1995; 
Petrosian  &  Liu  2004;  Liu  et  al.  2004b,  2006b)  and  implemented  in  the 
Stanford  unified  code. 


1566 


LIU,  PETROSIAN,  &  MARKKA 


Vol.  702 


produced  by  electron  heating,  could  contribute  to  generating 
seismic  waves  observed  in  some  flares  (Kosovichev  &  Zharkova 
1998;  Donea  &  Lindsey  2005).  (4)  In  the  long  run,  we  intend 
to  implement  time-dependent  particle  transport  calculation,  full 
loop  simulation  in  an  asymmetric  geometry,  return  currents,  and 
radiative  transfer. 

The  combined  treatment  of  the  particle  and  fluid  descriptions 
of  plasma  presented  here  opens  a  door  to  a  broad  range  of  appli¬ 
cations.  This  model,  with  proper  modifications,  can  be  applied  to 
environments  where  interrelated  particle  acceleration  and  trans¬ 
port,  and  plasma  flows  are  present,  such  as  (exo)planetary  au¬ 
roras  (e.g.,  Liu  &  Airapetian  2008)  and  flare,  on  other  stars  and 
in  accretion  disks  near  black  holes. 
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